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ABSTRACT 


An isoparametric finite element formulation solving a 
BSneral form of the transient field equation is presented in 
this thesis. The formulation is developed for fields in 
three-dimensional Euclidean space. A FORTRAN IV computer 
program employing double precision arithmetic, compact 
storage techniques, and providing the option of several 
numerical integration methods and types of finite elements 
is presented. Data input may be in rectangular, cylindrical, 
or spherical coordinates and variations in time integration 
step size are permitted. Time dependent boundary conditions 
are not considered. 

Comparisons of theoretical and computer solutions for a 
variety of test problems demonstrate close rohen. 


Instructions for use of the program are included. 
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IS SPDABESODUCTTON 


The 'quasi-harmonic' partial differential equation: 


9 дф 9 ὃ ϕ 9 [, 99 
Эх (x, 2). Зу (x, 23 + 52 ο >] 
+ [° =H ШЕ р = = 0 (1) 


ООШ ес to a variety of physical problems [Ref. 10]. Тһе 
Laplace equation is a typical and vell known example of a 
Specialized form of this equation. Since an analytical 
solution in closed form is seldom available for useful 
engineering problems, the development of a practical general 
purpose computerized solution method represents a useful 
addition to the analytical tools of the engineer. 

This thesis describes a general purpose FORTRAN IV 
computer program using the isoparametric finite element 
approach to solve the field eauation in general form. The 
program does not presently permit a time variation of boun- 
Bamyeconditions or a functional variation of material 
properties. 

The program allows selection of the linear, quadratic, 
or cubic isoparametric finite element with the option of 
specifying two through six Gauss points for numerical inte- 
gration in space. Data input may be a rectangular, cylin- 


drical, or spherical coordinates. Five time integration 
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routines are provided. Output is presented in a tabular 
form of node number, space coordinate, and functional value. 
The structure of the program is modular in order to 
accommodate improvements and adjustments due to new 
developments. 
This program was constructed and tested on an IBM-360/67 
computer system with OS/360 release 18. Other systems and 


installations will require modifications. 
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Ir. SEE ENT FORMULATION 


The 'quasi-harmonic' 


behaviour of some physical quantity, 0, is: 


9 9 Ó 9 
& s) x s 
2 
+ (0-0 80-р 26) о 


Commonly imposed boundary conditions are of the 


$ 7 j, 
39 οφ η 
κ. m T Е = s. = 8 К, 37 2 +: + op 


Where Pr represents the value of d specified on 


ft Eo and 2 are the direction cosines of the 


normal to the boundary surface; q is an outward 


unit area; Q a generation per unit volume; and Kur А КА; 


Q, K , K / k 


y, and p are material properties. 2 : 


may be functions of space and time. 


or field equation describing the 


(1) 


forms: 


(2) 


- 0 (3) 


the boundary; 
outward 
flux per 


2 


‚u, and» 


Reference 10 thoroughly develops the concept of a finite 


element and the process leading to a discretized form of 


equation (1) and the boundary conditions given by equations 


(C and (3). 


shape function, <N>, and nodal values of the variable, 
such that: 
Ó ‚ү {ф} 
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If the scalar variable, $, is described by a 


($), 


(4) 





the result 15: 
[H]{o} + [C] (0) + (6116) + {F} = 0 (5) 


[H], [C], and [G] are symmetric matrices defined, on the 


element level, as: 


Fel = m <N> πα ο ουσ (6) 
e 
[61° = f, «Ν᾽ ϱ «Ν» dx dy dz (7) 
e 
эм T ӘМ ON. „ӘМ 
[H] ° = f К <=> <— > + k <—>P* °? <— 
ν 9х 9 9 ΟΥ 
e y 
эм T ƏN 
. Κ, u^ a А6 У ασ (8) 
2 Ζ 


The load vector describing boundary conditions is: 


E a T T 
{F} = lis. Q<N> ` dxdydz + f. q <N> ds, 


e e 
+ Q <N>* α «ΝΣ as ) ШС (9) 
е 


where the τ are performed only over the surface or 
volume where the prescribed boundary condition applies. 
Element matrices and the load vector are assembled by 

standard methods to form equatron (5) which is thes cerc— 
tized form of equation (1). Adding the load vectors given 
by equation (9) to the rest of the nodal intensities 
guarantee the satisfaction of the boundary conditions given 
ENcuatron (3). Finally, the boundary condition of equa- 


tion (2) can be applied to the assembled set of equations. 
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o 


The finite element technique models a continuum as an 
Besembly of substructures; thus, for an element of small 
size, the parameters Kun Kye koe U, б, а, Q, ard q may be 
considered constant and the integrations required by equa- 


pens (6), (7), (8), and (9) reduced to an automated numeri- 


cal technique. 


TS OPARAMETRIC FINITE ELEMENTS 

The linear, quadratic, and cubic isoparametric shape 
functions are utilized in this development. The resulting 
elements contain 8, 20, and 32 nodal points respectively. 

A full development of the concept and properties of this 
family of finite elements is contained in Ref. 10; only the 
necessary results will be repeated here. 

A "local" non-dimensional system of reference (&,n,Zd) is 
used for convenience of calculation. Results are transformed 
into the "global" cartesian (x,y,z) system of reference prior 
to final assembly of problem matrices and vectors. The 
shape functions, М. (5,п,5), for the element types used in 
this development are given in Ref.-10 and need not be 
repeated here. The nodal numbering conventions employed for 
each element type are shown in Figures 1l, 2, and 3. 

Important results relating the shape function to the 


local and global reference systems [Ref. 10] are: 
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Figure 1. Linear Isoparametric Element. 
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Figure 2. Quadratic Isoparametric Element. 
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Cubic Isoparametric Element. 
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i 1 
οξ dx 
N, ΟΝ. 
In = [J] Jy (10) 
ΟΝ ΟΝ 
οζ 97 
where the Jacobian matrix, [J], is defined as 
οξ οξ 9 
- [9x  3y 92 
Ээх dy az 
οζ οζ οζ 
and 
dxdydz = det [J] d& dn dz (12) 


The partial derivatives of the shape function with 
respect to global coordinates are recovered from equation 


(10) by applying the inverse of the Jacobian matrix. 


B. MATRIX FORMATION 

Assuming constant material properties within an element 
wma lying equations (10), (11), and (12) to the matrix 
and vector formulations of equations (6), (7), (8) and (9) 
reduces the form of the integrals to: 


T 


νι det [J] dażdndac 


T 
y > < det (Il @eaednde +... . 


T 


Ie dnd 
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y Su бек [Т] ав 
S e 


ο. МЕ чес [у | 45 
© е 


where the integrations are performed in the local nondimen- 
sional coordinate system. The forms of equations (6), (7), 
(8), and (9) become ideally suited for a Gaussian numerical 


integration procedure. 


C. BOUNDARY CONDITIONS 

The specification of a value of the function 6 along a 
given boundary as in equation (2) is straightforward and 
requires no comment, 

Boundary specifications of equation (3) appear in the 
expression for a forcing vector given by equation (9). The 
first term of equation (9) containing Q represents a gener- 
ation rate per unit volume; the second term containing q 
represents a flux per unit area; and the third term contain- 
ing a represents a surface loss coefficient per unit area. 
Me terms containing q and Q contribute to a forcing vector, 
{F}, while the term containing a results in a correction to 
the corresponding element [H] matrix. 

If, in equation (5), q and a are zero and К, k „and к, 
are equal the boundary condition reduces to the requirement 
that the normal derivative be equal to zero. This condition 
is automatically imposed by equation (5) in the absence of 
any other boundary condition specification. 

Arpaci [Ref. 1] gives a number of boundary conditions 


applicable to the heat problem. Of these the Radiation, 
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Interface in Motion, and Change of Phase are not included in 
this development. 

Consider an element located at the surface of a body in 
ντε. Prom Figures 1, 2, and 3 observe that a corner node 
may be common to one, two, or three boundary surfaces per 
element and that a mid side node is common to one or two 
boundary surfaces per element. Each boundary surface may 
be subject to a different boundary condition specification. 
Tf two or more of the element surfaces have different boundary 
conditions imposed, a unique situation arises if the boundary 
Specifications include conditions on $ and q or а. In this 
case the boundary conditions must be specified in order. 
First, all boundary conditions involving q and a are formed 
and merged into the overall system of equations; then the 


boundary conditions involving $ are applied. 


D. SOLUTION TECHNIQUES 
For a fixed point in time and for steady state problems, 


equation (5) reduces to the form 
[Al{o} = {B} | (13) 


ООО (A) is an xn square matrix. The vector {B} is fully 
determined in all cases except where a nodal value of 6 is 
specified. If a nodal value, dir is specified the corres- 
ponding forcing vector term, Bi, will contain the unknown 


quantity of the equation. 
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M Stems of Linear Equations 


Equation (13) may be solved by standard methods of 
linear algebra. The method employed here to sort the known 
and unknown quantities of equation (13) was developed by 
Irons [Ref. 6] and has the advantages of preserving the 
symmetry of matrix [A] and permitting the determination of 
all unknown quantities of the system of equations. 

Meter applying the Irons [Ref. 6) modifications, 
solution of equation μία) is accomplished by the reduction 
of matrix [A] into the product of lower unit triangular, 
diagonal, and upper triangular matrices such that: 


EU — IL) [DJ [LJ 


Unknown quantities are found by performing a forward sub- 
STU lon followed by a backward substitution [Ref. 13]. 
Steady "State 
In steady state, equation (5) readily reduces to 
the form of equation (13) and is solved by the technique 
described previously. 
3. First Order in Time 
At fixed points in time, equation (5) reduces to the 
form of equation (13). The numerical integration schemes 
employed are straightforward and are UE below. The 
initial value of ($) is assumed known and time step size is 
denoted by h. 
a. Euler I Method 
This integration scheme is well documented by 


eh 


Crandall [Ref. 3]. The result, for the i point in time, 1s: 
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A WAN 


ο. een τε] (15a) 
mE Е соло (15b) 


b.  Trapezoidal Rule Method 
This approach to time integration is documented 


leet. 3 and Ref. ll. At mid interval 


(it = (or + to"? 
2 

(air da 10 - tor 
h 


When applied to equation (5), the resulting recursion 


formula is: 


1151 L 


1 qeptot*i e (tm + É [ο1){ϕ}᾽- 2{Ρ} (16) 


([H] + F7 


c. Linear Time Shape Function Method 
Reference 10 develops the following approach to 
time integration. Assume that {¢} may be considered as a 
product of the nodal values of {¢} at discrete times and a 


time shape function, M(t) , such that 
<M> = < > 
M M (t),M) (t) 

where, for a linear variation: 


|. (h-t) 
Mon h 


and 


M СВ 


1 


mae result is: 
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Applying this relationship and the Galerkin process described 


ih Ref. 10 results in the formula: 
2 1 s E E rA 
ИЕП P iC] = Е [Hi = p [C]) 19) ΤΙ; 


4. Second Order in Time 

At fixed points in time, equation (5) reduces to the 
form of equation (13). The integration methods employed are 
summarized below. Initial values of {¢} and are assumed 
to be known. Time step size is denoted by h. 

a. Euler I Method 

The Euler I integration method is well known and 

documented by many standard sources such as Ref. 3. The 


result, when applied to equation (5) is: 


оО τι ο στ.) (18a) 
en = "o + h ο... (18b) 
пог" (186) 


In the event that the [6] matrix equals zero these equations 
remain valid. 
b. Finite Difference Method 
Finite difference approximations are developed 
in Ref. 3. A four point backward difference approximation 


EOF (9) and a three point backward difference approximation 
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for {¢} have truncation errors of the same order of magni- 


tude, im and are: 


B e E (-(63373 & a(o)i72 - (9) 71420033). (194) 

eat? CTS (19b) 
Define 

{o*} = - Se? ο ας ο 509771, (196) 
and | 

{ф**} = = (5 πο cd (194) 


Equations (19a) and (19b) when applied to equation (5) give 


the following recursion formula 


2 3 ο. 
| 090 + 2 (El: 2 [СЈ) {Ф} = 
- {Е} - [Сб]{ф*) - [С]{ф**} (20) 


if [C] is not equal to zero and 


(8) € —- te) (27 » - (D - teltos) (21) 
h 
if [C] equals zero. 
Equation vas is known as the Houbolt method and is 
shown by Johnson [Ref. 7] to be unconditionally stable. 
Three starting values of ($) must be obtained by inde- 


pendent means such as the application of the Euler I method 


or a Taylor series expansion. 
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IEISEZEDESCRTFEITON OF THE PROGRAM 


Berinetional block ea of the program 1s presented 
in Figures 4 through 9. A modular structure was adopted in 
order to permit simple modifications to accommodate future 
improvements and the adaptation of the program for special 
purposes. 

The following paragraphs give a brief description of the 
function and operation of the program. Appendices A through 
D discuss the preparation of input data, list the important 
computer program nomenclature, and list the program. A full 
understanding of the operation of this program can only be 
achieved by a study of Appendices A through D and an under- 


standing of the principles and solution techniques employed. 


AR ASSUMPTIONS AND RESTRICTIONS 

The key assumptions and restrictions which impose 
restrictions in the type of problems that can be solved by 
the program are: 

l. Material properties Ko ES к, О, U, and a and 
problem parameters q and Q are constant functions of time 
and space throughout the element and are independent of the 
Blñetion é. 

2. Boundary values of ¢ are not functions of time. 


3. The principal axes of each element align witr the 


principal axes of the global coordinate system and are fixed 
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STEADY 


Figure 4. Subroutine Calling Sequence for 
Steady State Problem. 





MERGE (1) See Fig. 7 
| MERGE (2) See Fig. 8 


TIME11 
Or 
TIME13 






See Fig. 9 


Figure 5. Subroutine Calling Sequence for 
Pirst Order Time Problem. 
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| MERGE (1) [^ See Fig. 7 
| MERGE (2) | See Fig. 8 
| MERGE (3) Sed Pigi 8 


TIME 2 1 


OT 


——= See Fig. 9 


ENE 25 


NOTE: 


Figure 6. 





MERGEt2), 35 not called 1f the [el 
matrix is absent 


Subroutine Calling Sequence for 
Second Order Time Problem. 
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IMERGE (1) 


NOTE: STUFF not called for the Steady State Problem 


Figure 7. Assembly Routines for the Formation 
OF Matrix H]. 
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A рю 


MERGE (2) 
or 
MERGE (3) 






STUFF ¿olsa llaq 
from MERGE (2) for 
the Second Order 
Time Problem 






SHAPE | 


Figure 8. Assembly Routines for the Formation 
of Matrices [С] ала [С]. 


STEADY 
Tine | ου TIMES 


το ο or ЕМЕ 2 8 


BNSN | - ¢ MULT STUFFV 


NOTE: (1) STUFF called by TIME23 only 
(2) 4MUbT not called byssmsAnpDy 





Figure 9. Service Routines Used by 
Solution Routines. 
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in direction throughout the element. 


BE -imitations will be discussed in Section V. 


B. GENERAL PROGRAM ORGANIZATION 

The program may be divided into the three major indepen- 
dent functions: reading of input data, matrix formation, 
and actual problem solution. This organization is shown in 
Figures 4 through 9. 

All matrices involved in the problem are symmetrical. 

At the element level the complete matrix is formed; however, 
formation of the master matrices [H], [C], and [G] and prob- 
lem solution utilize a standard compressed storage mode to 
conserve computer core. The technique consists of placing 
the diagonal elements of each row of the matrix into the 
first column and retaining only the upper elements of the 
Marix. 

Reading of necessary input data is accomplished by sub- 
routine INPUT. Preparation of input data is discussed in 
Appendix B. An echo check of all input data is provided. 

Formation of the required matrices and vectors is 
controlled by subroutine MERGE. Matrices [H], [C], and [G] 
are formed sequentially. This permits output of the 
completed matrix to an external storage device and conserves 
core storage. Boundary condition calculations are performed 
during formation of the [H] matrix. Associated subroutines 


include CUBE, FORM, and JACOB. 
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The transient solution subroutines are similar in con- 
struction and differ only in the integration procedure 
employed. Time integration step size changes and data 
suut are controlled D» the solution subroutines. Certain 
functions common to all solution subroutines are performed 
by a set of service subroutines which include LDLT, BNSN, 


SEE. PUNC, and MULT. 


CS SUBROUTINES 
The MAIN program is simply a control segment that serves 
ES Call INPUT; MERGE (1), MERGE (2), or MERGE (3) in the proper 
order and as required; and the applicable solution subroutine 
ΠΥ, ΤΙΜΕΙΙ, TIME13, TIME21, or TIME23. No calculatıons 
are performed in the MAIN program 
e Subroutine INPUT 
This subroutine reads all input data and provides the 
appropriate echo check. The number of nodes per element is 
determined from the element type specification and, 1f 
required, cylindrical or spherical coordinates are trans- 
formed into rectangular coordinates by callina subroutine 
MARK. 
2. Subroutine MERGE (IMTRX) 
This subroutine forms the master matrix specified 
by the calling argument, IMTRX, by assembling each element 
matrix into the master matrix. A call to subroutine CUBE 
returns the desired element matrix.  IMTRX-1 corresponds to 


Ene TH] matrix, IMTRX-2 corresponds to the [C] matrix, and 
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IMTRX=3 corresponds to the [G} matrix. Subroutine BCOND is 
called during formation of the [H] matrix and, if required 
by the solution subroutine selected, the [H] or [C] matrices 
are placed into the appropriate storage location after 
formation. 
3. Subroutine BCOND 
This subroutine forms the element vectors and matrices 
required by the appropriate boundary condition specified and 
assembles the resulting vector or modification matrix into 
the forcing vector, {F}, or the corresponding [H] matrix. 
Formation at the element level is accomplished by a call to 
the appropriate entry point of subroutine CUBE. 
ae Subroutine STEADY 
This subroutine solves the steady state problem. 
After the application of appropriate nodal boundary condi- 
meee, a Call to subroutine LDLT and entry point SLV produces 
E j Lion. Data output is a printout of the (fé) апа (F) 
vectors and punched card output of the {¢ vector if requested. 
ЧООК Ὁκουιι TIME11 
This subroutine employs the Euler I time integration 
method, equations (15a) and (15b), to solve the first order 
time problem. One call to subroutine LDLT is required; sub- 
sequent calls to entry point SLV produce the desired solu- 
tion. Data output is a tabular presentation of the 
vector at specified points in time with punched card output 


of the {¢} vector provided af ee questee cr 


32 








E Subroutine TIME13 
This subroutine employs the "Linear Time Shape 
Function" as the time integration method (KINT=3), equation 
(17), or the Trapezoidal Rule time integration method 
(KINT=2), equation (16), to solve the first order time prob- 
lem. External core storage devices are utilized to store 
the [H] and [C] matrices in order to permit restarting the 
problem at each change in time integration step size. One 
call to subroutine LDLT is required for each change in time 
integration step size with the required solution obtained by 
sequential calls to entry point SLV. Data output is the 
same as for TIMEJL. 
TS byroutine TIME21 
This subroutine employs the Fuler I time integration 
EENM5d, equations (18a), (18b), and (18c), to solve the 
second order time problem. Operation of this subroutine is 
Identical to the operation of TIMEll except for the inclusion 
of the vector {¢} in the data output. 
8. Subroutine TIME23 
This subroutine employs a Four Point Backward Finite 
Difference (Houbolt) time integration method, equations (20) 
and (21), to solve the second order time problem. Starting 
values are obtained by a five term Taylor Series expansion. 
Operation of this subroutine is similar to the operation of 
subroutine TIME13 employing external core storage devices 
and saving the values of {¢} required to restart the problem 
after a change in time integration step size. Data output 


"e the same as for TIME21. 
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T. Subroutine CUBE(NPEL,NGP,ITYPE) 

This subroutine performs the Gaussian numerical inte- 
gration required for the formation of the various element 
matrices. Entry arguments are the number of nodes per 
element, number of Gauss points specified for the integration 
process, and the type of integration to be performed. Gauss 
Points and weighting factors were obtained from Ref. 9. A 
selection of two through six Gauss points with a default 
option of three Gauss points is provided. For each specified 
Gaussian integration point, the numerical evaluation of the 
integrand is determined by calling subroutine FORMC at the 
appropriate entry point and applying the proper weight 
factor to the result prior to summation. 

The specialized volume and surface integrals required 
for boundary condition calculations are obtained by entering 
P subroutine at entry point CUBEB. 

10. Subroutine FORMC (XI,ETA,ZETA,NPEL) 

This subroutine performs the products of «му! >, 
«5 ie eu >, etc. required at each point in local coor- 
dinate space for the numerical integration process. Values 
of the shape function and partial derivatives of the shape 
function with respect to the local coordinate system are 
determined by a call to subroutine SHAPE. Calling subroutine 
JACOB returns the required values of [J], [J] апа Det [J] . 
Necessary entry arguments are the local coordinates of the 
point in space under consideration and the number of points 


per element. Various entry points are employed to accommodate 
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the specialized calculations required in the formation of 
each element matrix or vector and volume or area integrations. 
Πα. =sübroutine SHAPE(XI,ETA,ZETA,NPEL) 

This subroutine forms the shape function and the 
partial derivative of the shape function with respect to the 
local coordinate system at the point in local coordinate 
| space specified by the entry argument. The type of element 
to be considered is determined by the specification of the 
number of points per element in the calling argument. 

IE Subroutine JACOB(NPEL,AJ,DCOL,CORD,AJI,DTJ) 

This subroutine evaluates the Jacobian matrix, 
determines the inverse of the Jacobian matrix and the deter- 
minant of the Jacobian. The number of points per element, 
global coordinates of the node points, and shape function 
partial derivatives with respect to the local system of 
reference are required entry arguments. [J] is returned in 


a (a): 


In AJI, end det[J] in BJ. 

The specialized calculations required for a surface 
integration are accomplished by entering this subroutine at 
Eutrv pornt JACOBB. | 

E Subroutine LDLT @,N,M,ASN) 

This subroutine reduces the matrix [A] into the 
products of lower unit triangular, diagonal, and upper tri- 
angular matrices, [L] [D] pude. , required by the process of 
solving the linear system of equations described by 


[Al {x} = {B}. N and M are the number of equations and half 


band width respectively. ASN is an arbitrarily small number; 
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numbers less than ASN squared are considered to be zero. 
ο τν point SLV(Z,B,N,M) performs the forward substitution 
and backward substitution required to solve the system of 
equations. 

14. Subroutine MARK (C,K,M) 

This subroutine converts coordinates given in a 
cylindrical or spherical system of reference to a rectangular 
System of reference. C is a matrix containing the specified 
coordinates, K is a code to indicate the system of reference 
employed, and M is the number of coordinates specified. 

u Subroutine BNSN(A,N,APZRO,ABIGN) 

This subroutine determines the arbitrarily small 
number required by subroutine LDLT and the arbitrarily large 
number required by the equation solving process. APZRO is 
determined from the matrix [A] by finding the average value 


Me he diagonal elements and dividing the result by 101°. 


ABIGN is m x ASN. This provides the numerical spread 





necessary for the IBM/360 computer installation emploved for 
{this program development. 
Ecc Subroutine MULT(A,B,C,N 
This subroutine performs the matrix multiplication 
[A] 7 = {C} in a compressed storage form. This multipli- 
cation is required by the various time integration subroutines. 
N and M are the number of equations and half band width 


respectively. 
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ie Subroutine STUFF (A,B,N,M) 

This subroutine places the matrix [A] into matrix 
[B]; or entry STUFFV(AA,BB,N,MC) places the vector {B} into 
the specified column, MC, of matrix [AA]. 


ШООК oubroutine PUNC(A,B,N,K) 
This subroutine produces the punched output requested 


SE the Solution subroutines STEADY, TIME11, TIME13, TIME21, 


| BE TIME23. 
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IV? TEST PROBLEMS 


The test problems employed here are designed to test and 
verify the operation of the computer program. The four cate- 
gories of tests performed were assembly, steady state, first 
order in time, and second order in time problems. 

The assembly tests served to verify the calculations and 
integrations required to form the problem matrices and vectors 
on the element level. Steady state problems were used to 
verify the matrix assembly process, linear equation solution 
process, and the validity of boundary condition formulations. 
Transient problems of first and second order in time served 


ΠΝ... τ the various solution subroutines utilized. 


"ῬΜΕΤ,Υ TESTS 

A rectangular region in space provides a convenient means 
of verifying element matrix and vector assembly procedures. 
The calculation of all problem matrices and vectors is 
stralghtforward and element calculations of the [H], [C], and 
[G] matrices were verified in this manner. Similarly, the 
basic assembly processes of subroutines MERGE and BCOND are 
readily verified. 

The integrations required by equations (6), (7), (8) and 
(9) are performed by a Gaussian numerical integration scheme. 
In the Gauss procedure, n sampling points will exactly 
integrate a polynomial of degree (2n-1). The integrations 


necessary may involve polynomials of degree 2 through 14 in 
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SI En direction of the local coordinate system. Poly- 
iomial degree will depend on the shape of the element 
GEE. For practical applications, a selection of two 
ΠΝ. πχ Gauss points is sufficient. For a rectangular 
hape the linear quadratic, and cubic elements are exactly 
ntegrated by 2, 3, and 4 Gauss points respectively.  Veri- 
'ication of the number of Gauss points required to exactly 


ntegrate elements of other shapes was not considered. 


ZEN STATE PROBLEMS 
NAL Temperature Distribution in a Cylinder 
Consider a cylinder of the geometry specified in 
"Figure 10. If the sides are perfectly insulated a linear 
emperature distribution results when end conditions are 
feted [Refs. 1 and 5]. Specify the thermal conductivity, 
20" BTU/nx/ft/*F and divide the region into two finite 
>lements. The top surface is subjected to an anbient tem- 
mature of 100°F. 
If the bottom of the cylinder is subjected to a 
laut flux, aq, of 500 BTU/hr/ft* or maintained at 300°F with 
Sat transfer coefficient, a, of 5 BTU/hr/ft^/9F identical 
cemperature distributions result. Theoretical and finite 


2lement solutions are shown in Table I. 
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Figure 10. Axial Cylinder Problem Geometry. 


TABLE I 


AXIAL TEMPERATURE DISTRIBUTION IN A CYLINDER 





Temperature (°F) 


ا ر 


Distance 2 Quadratic 2 Cubic 
from Theoretical Element Element 
Bottom(ft) Solution todel Model 
0.00 200.00 200.00 200.00 
ЧОО? PESAS ο 98 
1660 0 175.00 175.00 
E33 166.67 που οἱ 
2800 15000 150. 00 150.00 
2.67 IS en. 33 
500 25400 125.00 
6233 DEG O7 ΠΗΡΕ 
4.00 100.00 100.00 100.00 


e Y n‏ 55 ا 
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EE пуат Temperature Distribution in a Cylinder 


Consider a hollow evilinder “ofermii nite length” A 
Begment of such a cylinder is shown in Figure ll. Heat is 
ost by convection through the outer vall. Consider the case 
f the interior surface maintained at a constant temperature 
r the case of a specified heat generation within the cylinder 
With the interior surface perfectly insulated. Assumed 
roblem parameters are: 

a. Inside radius (ri) one 

b. Outside radius (τω) 9 7. 


C. Heat transfer coefficient (a) 36 BTU/hr/ft^/9*F 


ο Thermal conductivity (K) ΘΙ“; πε; E 
e. Ambient temperature OT 

f. Inside wall temperature σος TE 

g. Heat generation rate (0) 22 BTU/hr/ft? 


Theoretical solutions to this problem are readily 
obtained by methods outlined in Refs. l and 5. Comparisons 
Of solutions obtained by finite element models and theoreti- 
cal means are given in Tables III and IV. 

A finite element model of nr problem is constructed 
by assembling elements radially to form a segment of the 
cylinder. Since the problem is symmetric and no heat is 
conducted perpendicular to any radius, any convenient 
specification of element arc width, 9, and axial thickness, 
L, may be used. A small element arc is desirable since the 


isoparametric finite element used will provide a better 


approximation to a true circular arc. A large arc width 
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will introduce some distortion in the temperature distri- 
bution. Note that the sum of the forcing vector elements 
kepresents the total heat flux through the body in this case. 


The various finite element models used are given in Table II. 


TABLE II 


FINITE ELEMENT MODELS FOR RADIAL CYLINDER PROBLEM 


Model 

D L-24 0-6 COMI 
Element 
Thickness (L,in) 1 1 1.5 1.5 
Element 
Width (9,2) 2 2 6 6 
m Linear Linear Quadratic Cubic 
Type 
Number of 
Elements 12 24 6 4 
Number of 
Nodes 32 100 80 92 
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Figure 11. Radial Cylinder Problem Geometry. 


PABLE III 
RADIAL TEMPFRATURE DISTRIBUTION IN A HOLLOW CYLINDER 
КИП THE INTERIOR SURFACE MAINTAINED AT 1500 °F 
Temperature (°F) 


Radius Theoretical Model Model Model Model 


(in) τε L-24 Ab Cra 

3 1500.00 159000 1209200 1500. 00 1500200 
¿ 15898. 54 1188585 55 188.57 153 
5 946.96 947.34 9717.06 947.01 946.99 
6 74956 749.94 55 749.63 749262 
7 382.67 583300 502577 οσο ο 582.74 
8 438.11 430.38 438. 20 438. 19 498,19 
9 220,99 310. 51 3102767 310769 3102008 


Heat Flux (q) 
Theoretical 9.45 9.45 85.03 85.03 
calculated 9.46 9.46 85.03 85.02 
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TABLE IV 


RADIAL TEMPERATURE DISTRIBUTION 
WITH INTERNAL HEAT GENERATION 


Radius Theoretical 

(in) 

3 175,28 
4 172.00 
5 1662777 
6 η 
И 144.35 
8 128.60 
9 1O00 


Heat Flux (q) 
Theoretical 
Calculated 


IN 


A HOLLOW CYLINDER 


Temperature (°F) 


Model 


El 


L7 
17954 
166. 
151. 
144. 
125. 


FO, 
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2 


34 
03 
ὦ: 
3 
34 
60 


99 


291 


97 


Model 


Q- 


ΠΡ 


1735 


166. 


144. 
1,28. 


110. 


14. 


14. 


6 


30 


02 


[9 


2 


36 


62 


02 


14 


14 


Model 
C-4 


2732 
1738 
166. 


15 


144 


шав 


1.1708 


148 


14. 


30 


02 


79 


14 


30 


62 


02 


14 
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PR Chimney Problem 


Consider the two dimensional heat transfer problem 
posed by a rectangular chimney with the inside surface 


Mam ained at lOO°F and the outside surface maintained at 


ОЕ. The steady state temperature distribution may be 


obtained by constructing a finite element model consisting 


of the desired number and type of elements assembled to 
form a plate having the required geometry. The geometry of 
Buea model utilizing symmetry boundary conditions is shown 


in Figure 12 and the results of various finite element models 


are compared with a solution obtained by a relaxation tech- 


nique in Table V at the points shown in Figure 12. 
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Figure 12. Rectangular Chimney Problem Geometry. 


TABLE V 


TEMPERATURF DISTRIBUTION IN A RECTANGULAR CHIMNEY 





Temperature (°F) 





5 Cubic 10 Quadratic 40 Linear 





Ent (x,y) Relaxation Element Element Element 
Method Model Model Model 
a з 257.5 44.3 8572 
b 2 368:9 due 5 35 3 
E Il 2026 194 203 
d 251 39/2 aged ο. ο 38.7 
e Salt 47.7 47.0 476 5 
τ 4,1 49.3 48.7 48.9 49.1 
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err Rol ORDER TIME PROBLEMS 

Transient heat conduction in various bodies is a common 
form of the field equation that is easily modeled by the 
finite element process. Transient heat conduction ina 
hollow cylinder and a prolate spheroid will be considered 
Here. 

Euler I, Trapezoidal, and Linear Time Shape Function time 
integration methods were employed. Each method yielded 
excellent results vhen used with an appropriate time step 
size. However, with the cylinder problem, the Euler I 
method was unstable for large time increments. The Trape- 
zoidal and Linear Time Shape Function methods produced 
excellent results in all cases. 


l. Transient Temperature Distribution in a Hollov 


e tinder of Infinite Length 


Consider the hollow cylinder of Steady State Test 2, 
Figure 11. Thermal diffusivity, a, is specified as 
E25 ИО hr where a = k/pc. Initially the cylinder is at an 
ambient temperature of 70°F; suddenly the inside wall tem- 
perature is raised to 1500°F and maintained at this temper- 
ature. All other conditions are the same as for the Steady 
state Problem. 

An approximate analytical solution, which is valid 
until the temperature wave propagates to the outer wall of 
the cylinder, can be found by applying methods of solution 


ENScussed in Ref. 1. The result is: 
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r-r. 
= (To η vr,/r exte [ 3r i? 
Dre 


The finite element models employed here are identi- 
cal to the models employed for the Steady State case. Tran- 
Sient solutions obtained by the Trapezoidal and Tine Shape 
Function methods and using the linear element model converged 
to the same results obtained in the Steady State Problem for 
corresponding model. Solutions with Quadratic or Cubic 
models and Euler I integration were not pursued to steady 
state. 

bles VI through IX compare analytical and finite 
element solutions at various problem times for the 100 node, 
24 Linear element model used with the Steady State Problem. 
Comparable results were obtained with Quadratic and Cubic 
element models. The following abbreviations apply to the 


time integration methods: 


ου E-I Puder I 
b. TRAP Trapezoidal 


En LTS Linear Time Shape Function 
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TABLE IX 


RADIAE TEMPERATURE DISTRIBUTION 
CYLINDER AT 1.50 MINUTES 


Temperature (°F) 


Method 
Bact B-l TRAP 
h (min) „000 5m w00 pa 
Radius (in) 
3.00 roo. 00 rome. 00 17500: 00 
bu25 1240.81 1245.47 1246.07 
325850 TUNE TIS lI0T6.54 1077.44 
JO 506291 61506 846.53 
4.00 634.79 642.85 644.53 
5.00 2 0569 224517 225454 
6.00 5» 95297 96.47 
7.00 2869 ο ο) 7286506 
8.00 70.17 70.14 O15 
9.00 7001 /0 01 O01 





he 
SUS CO 1 min, .001 thereafter 
ХЖ 
mee to l min., .0l thereafter 


I2 


LTS 


ΟΦ 2u: 


1500. 
1246. 
ТОЧ 
816. 
644. 
2:254 
gs 
728 


0. 


00 


00 


33 


38 


38 


50 


50 


68 


15 


7001 


IN A 


TRAP 


zu 


1390 200 
1246.00 
1017. 32 
E6. 37 
644. 35 
22539 


96.42 


70425 


7001 





LTS 


ol 


1500. 
1245. 
1017. 

816. 


643. 


ος. 
2 
70. 


70. 


00 


84 


04 


01 


95 


29 


50 


70 


185 


01 





fee transient Heat Conduction in a Prolate 


Spheroidal Solid 


An analytical solution for transient heat conduction 
in a prolate spheroidal solid is developed in Ref. 4. 
Results of a finite element solution to this problem are 
presented in Refs. 10 and ll. 

A non-dimensional finite element model using one 
cubic element or four Quadratic elements was constructed. 
The Linear Time Shape Function of time integration was 
employed giving results comparable to the results of Refs. 
10 and 11. A comparison of theoretical and finite element 


results is presented in Figure 13. 


NN SECOND ORDER TIME PROBLEMS 

The transient vibration of a perfectly elastic membrane 
1S a two-dimensional problem that can bo formulated as a 
finite element model of the field equation. 

Consider the square elastic membrane of length a and 
specified material properties Tor Pr and v. p is the 
Enitial tension of the membrane, p is the surface density, 
V is the damping factor, and СА = T. /p. The partial dif- 


ferential equation describing the motion is: 


2 2 2 
Dg u 2 0U e 


9х ду t 


а? 
f 
а? 
с 


C 


ره 

rt 
ره‎ 

NO 


where u(x,y,t) is the displacement of the membrane. 
An exact mathematical solution is derived in Ref. 2 and 
may be specialized to a simple closed form if an initial 


displacement of the form 
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i: les im 
u(x,y,0) = d sin = sin > 


and boundary conditions of simply supported edges are assumed. 
Two finite clement models were constructed. One models 
the entire membrane as a 36 Linear element plate; the other 
as a 9 Quadratic element quarter plate imposing symmetry 
boundary conditions at the centerlines. Euler I and Finite 
Difference solution techniques were employed with a time 


step size of .001. Assume the following material properties: 


a. a= 1.5 
CTE 
ce VY = 3.0 
BS - .25 


Both time integration methods produced good results; how- 
ever, the Euler I integration method was unstable at large 
time increments. 

1. Undamped Membrane 


The analytical solution is 


15 à TX ° TY 
Se πμ μα + 
a a 





u = d cos 


Results obtained with the finite element models are compared 
to the analytical result for the center node displacement in 


Table X. 


D. Damped Membrane 


The analytical solution is 


-vt 


νη. TX TV 
u = d e (cos wt + P sin 0t) Sin = “os 


a 
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CENTER NODE DISPLACEMENT OF AN UNDAMPED MEMBRANE 


TABLE X 


Displacement 


EULER I Method 


ЕНЕ ТЕ ΡΙΞΤΕ 





Method 

Time Pret Linear Quadratic Linear Quadratic 

Model Model Model Model 
.00 "2500 .2500 .2500 .2500 .2500 
05 2378 LIZ ο . 2370 2868 
.10 ‚2023 .2007 .2108 . 2002 .2014 
215 . 1469 . 1439 ‚1456 ‚1435 . 1453 
. 20 10773 .0726 20763 .0723 „Оо 
E25 .0000 .0061 = 0020 οσο =. 0018 
"30 20773 .0844 — ¿QUIEN 0889 -.0784 
E35 po MIO .1543 —. 1496 189552 -.1483 
. 40 = 2023 ZOO -.2041 . 2073 “ΜΘ 
«45 5 .2428 = 12—07 ‚2406 шы оэ 
850 PASO 523 “ΕΠ . 2498 -.2495 

where 
ш = N: - у? апа Ne = т/а 


Results obtained with the finite element models are compared 
to the analytical result for the center node displacement in 


Table XI. 
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TABLE XI 


CENTER NODE DISPLACEMENT OF A DAMPED MEMBRANE 


i nn nn 


Displacement 


ne nn —_——— 


` J 7 : STONI 4 
FULER 1 Method FINITE DIFFERENCE 


Method 

EE . Zins eL кз оо 
к 0 Qo ——— — V n n 
EDO . 2500 .2500 22200 .2500 22300 
E05 82389 22:94 s2382 ¿2982 2300 
‚10 22105 20) «2100 «2086 2097 
1.5 .1716 411092 12.05 “БО. „1704 

. 20 283 1249 ZO +1250 . 1274 
"25 .0854 .0813 „0841 .0815 .0843 
30 .0464 .0418 .0454 .0422 .0457 
5 E0135 „0089 -0122 -0095 ο.» 
.40 2:01:19 = 0 02 =, O12 8 -.0156 - OR 
.45 0298 ΘΒ = 0609 = 0326 πα αμ 
850 -.0405 -.0434 -.0411 EEE -.0406* 
55 -.0450 m =. 0558 20520 
.60 -.0445 | πη MIGO) 
"65 -.0404 -.0621 -.0644 
2570 -.0340 ΠΡ OSD 
m5 -.0263 0457 _.0490ї 
— | |] | | | ||. 1] Qo LS SL LLL —-— 
* 

Changed step size to .0l 

Ке сеа si size ton. 01 at time = xb; Solution 


pecame unstable 
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V. CONCLUSIONS AND RECOMMEPDATIONS 


The computer program presented in this thesis provides 
an accurate and reliable means for solving a variety of 
problems. The use of this program and efforts to increase 
its versatility are highly encouraged. 

The Linear Time Shape Function and Four Point Backward 
Finite Difference (Houbolt) time integration methods provide 
the most consistent ane reliable results. Investigation of 


other time integration schemes should be pursued. 


A. RECOMMENDED MODIFICATIONS TO THE PROGRAM 

A number of modifications to the present program are 
desirable in order to increase the versatility £ the program. 
Most of these modifications are straightforward and are the 
result of the simplifying assumptions made during develop- 
ment of the program. The recommended modifications should 
Bresent no major programming problems. 

1. Anisotropic Materials 

Расіі су for coordinate axis rotation will provide 

the program with the ability to accommodate problems in 
which principal material axes do not align. This may be 
accomplished by the addition of an appropriate subroutine 
to perform the axis rotation. This subroutine should be 
called from MERGE at the time of element formation. Minor 


Changes to subroutine INPUT will be required. 
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Time Dependent Boundary Conditions 


mie ability to account for a time variation in the 
boundary value specified for {¢} is highly desirable. In 
meanmciple the modifications required are straightforward. 
Mee required modifications will only affect the solution 
u os TIMBll, TIME13, TIME21, and TIMF23 if data input 
is accomplished by providing updated values of {¢} and {6}, 
if required, in tabular form. The modifications will require 
attention to detail and a thorough understanding of the Irons 
technique as applied in the program. At present, boundary 
values of a) are considered to be zero; this effect on the 
solution matrix corrections and the effect of time dependence 
πο (F) vector must be carefully considered. 

Time dependence of material properties implies the 
recalculation of problem matrices at each station in time 
and possible iteration. Considerable computational oo 
is required and efficient programming techniques are 
necessary. 

3. Non Linear Material Properties and Boundary Conditions 

A functional dependence of material properties and 
non linear boundary conditions such as the radiation boundary 
ИЯ oí the heat problem are real physical situations. 
Mme ability of the program to accommodate this situation is 
needed and feasible by various iteration schemes. 

4. Problem Size 
The program is presently designed for problems of up 


to 110 equations with a maximum band width of 64. Practical 
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E ical problems will require much more computer storage 
E. Modifications to accommodate more equations in core 
require a careful adjustment of the appropriate common and 
dimension statements and attention to the equivalence state- 
ments in the time integration subroutines. 
specialized Data Input and Output 

Hand preparation of input data is tedious and subject 
to numerous errors. For practical physical problems of even 
moderate size, automated data input is a necessity. 

Each class of physical problem will require data 
ωρα in a specific form. Considerable attention should 
be given to the development of various subroutines to 


tailor data output to a given class of problem. 


BEZZRECOMMENDED FUTURE STUDIES 

This thesis was concerned with the development of a 
Computer program with the capability of solving a class of 
physical problems. This objective was accomplished and 
recommendations for improvements stated. Questions raised 
during the course of this development and beyond the scope 
of this work are: 

1. Convergence criteria of the various time integration 
5, particularly the Euler I technique. 

2. Selection criteria for finite element type and, in 
particular, the type of element best suited to a given 
EESuetry or physical situation. 

3. Adaptability and usefulness of other finite element 


types. 
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APPENDIX A 


COMPUTER PROGRAM NOMENCLATURE 


A complete listing and description of all variables used 


metic program is not practical. The items listed in this 


appendix are common to several areas of the program and will 


assist the reader in a study of the program. For convenience 


items are listed by variable type. The definition or func- 


MON and dimension, if applicable, of each item is given. 


А. INTEGER CONSTANTS 


IGO 


IPUN 


IP ROB 


HERT 


KBS 


KBT 


KC 


KINT 


MRH 


NBAND 


NEL 


NELBC 


NGP 


NMAT 


NNBC 


Element type 
Punched output request 
Problem type 


Number of time increment steps to skip 
between print out of results 


Boundary condition surface 

Boundary condition type 

Coordinate system of reference 

Time integration type 

Time step size index 

System band width 

Number of elements 

Number of element boundary conditions imposed 
Number of Gauss points 

Number of materials 


Number of nodal boundary values specified 
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NNODE 


NPEL 


NTSC 


NTP 


Number of nodes 
Number of nodes per elerent 
Number of time step size increments specified 


Number of times {¢} is to be printed 


КО REAL CONSTANTS 


ABI GN 

AKX 

AKY 

AKZ 

APZRO 

DT 

DTJ 

EACT 

EHIREF 

BEIM 

WFACT 
WX,WY,WZ 
EL, ETA,ZETA 
11,21 or 
Геб ҮҮ, 272 

EE VECTORS 
COL (32) 
ENDT (5) 
DTIM (5) 

TE (32) 


110) 


Arbitrarily large number 

[HJ matrix property in the x direction 
[HJ matrix property in the y direction 
[II matrix property in the z direction 
Calculated value of zero 

Time step size 

Determinant of the Jacobian 

[elmer [G] matrix material property 
Reference value of % 

Initial problem time 

Products of Gauss weight factors 

Gauss weight factors for specified direction 
Local cooraanates t, n, amd c 


Global coordinates x,y, and 2 


Element level storage of the shape function 
End time of the various time integration steps 
Time increment for each time step size 

Forcing vector at element level 


Problem forcing vector 
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(5) 
KBCOND (48) 


KEL (48) 


KNODE (110) 
KSURF (48) 
PHI (110) 
PHIDOT (110) 
BHTDDT (110) 
PHIWRK (110) 


TIME (6) 


E MATRICES 


E , 3) 

Eur (3,3) 
BGC (110,64) 
BGG(110,64) 
BGM (110,64) 
COARD (110, 3) 
CORD (32,3) 


EOL (32,3) 


ο 52,32) 


ENW(32,32) 


NCON (48,33) 


PROP (10,5) 


OUT (110,6) 


IPRT for each time integration Step size 
KBT for each element 


Element number for application of boundary 
conditions 


Node number for specified ($4) 
KBS for each element 

(ó) vector 

(6) vector 

($) vector 

Working vector 


Time label for output 


[J] matrix 


E 


[J] matrix 


[C] matrix 

[G] matrix 

[C] or [H] matrix, depending on problem type 
Global coordinates 

Global coordinates for element calculations 


Derivatives of the shape functnon with respect 
to local coordinates. Column I is 3/35; 
Column II is 3/ə3n; Column III is 9/9r. 

Work matrix for accumulation af [H] , es s 


er [G]; 


, 


Work matrix for element level calculations 


Connectivity matrix for each element, column 
1 1s material number 


Haterial properties 


Equivalenced matrix for storaqe of output data 
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Gauss points 


Gauss weighting factors 
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APPENDIX B 


MIS PRUCTIONS FOR PROBLEM DECK PREPARATION 


Although the procedure is straightforward, preparation 
OE input data for the program requires meticulous attention 
Fo detail. Errors are easy to make and difficult to locate. 

iat data preparation should, in general, follow the 
steps outlined below before attempting to punch data cards. 
ISS Of the standard Fortran Eighty Column Coding Sheet 
is highly recommended. A sample problem deck follows this 
appendix. Integer constants must be right adjusted in the 
appropriate field. 

Preliminary preparation: 

1. Formulate the problem in terms of equation (1) and 
boundary conditions in terms of equations (2) and (3). 
Identify the appropriate material properties. In the case 
of a one or two-dimensional problem, certain properties 
must be specified as 1.0; a specification of 0.0 will lead 
to a singular system of ae 

2. Divide the continuum into the desired number and 
type of elements. Observe the restrictions of 48 elements, 
110 nodes, 64 bandwidth, 10 materials, and 48 element level 
boundary conditions. 

3. Assign node numbers as Васа, observing the band 


math restriction. 
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4. List the element connectivity in accordance with 
EC conventions of Figures 1, 2, or 3. View £, п, апа б 
AS x, y, and z respectively. Any other convention will 
lead to serious error. 

oe List all nodal points and coordinates in any con- 
venient order. 

6. List the applicable boundary conditions and the 
appropriate codes and material properties. 

mee τί applicable, select the solution subroutine 
desired and list the required time integration parameters. 
Mor the Houbolt method the specification of ENDT must allow 
enough time to generate the starting values required for 
the next step in time. 

Preliminary preparations are now complete; precise card 
Bonching instructions follow. Card format is given in paren- 
thesis followed by specific instructions where necessary. 
Recall that blank columns are read as zero. 

Be Card (10А8) 

Eighty columns on one card for any desired title 
NEO rmation. Column 1 of this card will not be printed; 
if a l is punched, the output starts on a new page. 

2. Problem Parameters (1015,10X,F20.0) 

One card containing the following data: 
a. Cols 1-5 : Number of nodes (NNODE) 
b. Cols 6-10 : Number of elements (NEL) 
c. Cols 11-15 : Number of materials (NMAT) 


d. Cols 16-20 : Number of Gauss points (NGP) 
Specify 2 through 6; default to 3 
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ας, 


k. 


Cols 21-25 : Element type 


l^ Linear element 
2 = Quadratic element 
3°, Cube "element 


Cols 26-30 : Number of element level 
boundary conditions (NELBC) 


Cols 31-35 : Wümber of nodal boundary 
values specified: (NNBC) 


Code 36-40 : Próblem type (IPROB) 


l = Steady state 


2 = Eirst order in «Mane 
3= Second order in time 
4 = Second order in tie without {С) mataran 


Cols 41-45: Coordinate system of reference 
(KC) 


1 = Rectangular coordinates 
2 = Cylindrical coordinates 
3 = Spherical coordinates 


Cols 46-50: Punched output requested (IPUN) 
0 = No punched output desired 


l Punched output desired 


II 


Cols 61-80: Reference value of 6 (PHIREF) 


3. Material Properties; One card per material for the 


Steady state problem or two cards per material for the 


transient problem. 


a. 


Card I Jug 3E20.0) 
Cols 1-10: Material number 
Cols 11-30: k 
x 
cols 31=50: κ 
y 


OM 





ο: ο) τα К, 
b. Card II (110,2F20.0) 
Cols 1-10: Material number 
Cols 11-30: [G] matrix coefficient 
Cots 31550: 9 πο τὸ coeftfrcyent 
' 4. Node cards; One set of cards for the steady state 
Problem or two sets of cards for the transient problem. 
Each set has NNODE cards. 
a. Sw I (IIO,3F20.0) 
Cols 1-10: Node number 


Cols 11-30: x coordinate or radius 


Cols 31-50: y coordinate or angle in degrees 
relative to the x axis (+CCW) 


(ls 51-70: æ coordinate or angle in degrees 
relative to the z axis (+CW) 


pe Seb rr Chi 289050) 
Cols 1-10: Node number 
Gols vi l-30 re inieral value or m 
Coles] -SOsinrerat varie or Ó; 
Leave blank if first order time problem. 
5. Node boundary values; NNBC cards required but omit 
E IMEC = O (110,F20.0) 
Cols 1-10: Node number 
σοις 11-30: Specified value of Ф: 
5. Connectivity data; NEL cards required for the Linear 
Or Quadratic element and NEL x 2 cards required for the Cubic 
element (215,32I3) 


as Card 00275 ,2683 ) 
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Cols 1-5: Element number 

Cols 6-10: Material identification number 

Cols 11-79: Connectivity for nodes 1 through 23 
b. Card II (913): Use with Cubic element only 


Cols 1-27: Commectiviry «for cubic element 
nodes 24-32 


ΙΙ... Boundary condition data; NELBC cards required but 
meet if NELBC = O0 (315,F20.0) 
Cols 1-5: Element number 


Cols 6-10: Boundary condition type (KBT) 


] = Volume integral of loading (Q) 
2 = Surface integral of loading (q) 
3 = Surface loss coefficient (a) 


Cole 11-15: Wpe of #ntegral and surface 
applicable (KBS) 


O = Volume integral; Use with KBT=1 only 
11 = Surface integral over &-* ] face 
+2 = Surface integral over n=+ 1] face 
+3 = Surface integral over I=+ ] face 


Note: KBS = 1,2, or 3 must be signed plus (+) 
or minus (-) to indicate the 
applicable integration surface 

Cols 16-358 Boun@@ry comdition prope@ty value 

8. Transient problem data; NTSC + 1 cards required but 
omit for the steady state problem. 


um Сата т, опе card (2710. 0ος Ε20.0) 


Cols 1-10: Number of changes in time step 
size (NTSC) 


Cols 11-20: Time integration method (KINT) 
Default to method 3 
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1 = Euler I method 
2 = Trapezoidal rule method 


3 = Linear Time Shape Function or 
Backward Finite Difference method 


Cols 31-50 area) problem time 
Im Set. ll, NISC cards (Т10,2Е20.0) 


Cols 1-10: Numaer of time intervals to skip 
before printing output (IPRT) 


Cols 11-30: Time interval (DT) 


Cols 31-50: Tote] problem time to ema of 
interval (ENDT) 


Data preparation is now complete. Cards should be 
arranged in order of ascending node numbers within each set 
er group and each set or group of cards stacked in the order 
Brepäred. Refer to the sample problem deck following this 
appendix. 

Several problems may be stacked and submitted if desired. 
In order to terminate a given computer run add one last 
Mata card with the word "STOP" punched in columns 1-4. 
Standard JCL cards and deck structure is shown in Appendix 
C. Print out of the program is suppressed by the "//FORT. 
SPRINT DD DUMMY" card; if a program listing is desired 
remove this card. The JCL card requesting additional output 
should be removed unless voluminous output is expected. 


Good Luck! 
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